Muscle_RNAseq
Sam
11/02/2022
Muscle RNAseq results
Methods
Differential gene expression analysis was conducted using a multilevel or repeated measures design. The repeated design refers to the fact that multiple measures were performed in the same individuals, so we have a paired comparison experiment for each genotype group allowing us to compare the effect of treatment, in this case euglycaemic clamp, between the same individuals. However, additionally we would like to compare the effect of treatment between genotype groups. This was perform using DEseq2 a model matrix was constructed using the following formula: ~Genotype + Genotype:individual number + Genotype:Stimulation. Comparsions between groups were performed at basal level using controlling for age + sex + BMI. Functional enrichment of genes was performed using gprofiler. Genes were ranked based on the absolute LogFC and submitted as an ordered query giving more weight to those genes with higher absolute logFC. LogFC changes we shrunk using the apeglm package.
Number of significantly differentially expressed genes
## 1651 siginificant genes in WT between Insulin & Basal## 917 siginificant genes in MUT between Insulin & Basal## 0 siginificant genes between Basal WT & Basal WT## 0 siginificant genes between Insulin MUT & Insulin MUT
Overlap of significant genes with old model non-repeated measurements design
## Correlation of significant genes between old and new model 0.9584546WNT genes
Genes related to WNT signaling are differentially expressed, most notably, BCL9, BCL9L, WNT11, WNT9A, WNT4.
BCL9, together with its paralogue gene BCL9L (BCL9 like or BCL9.2), have been extensively studied for their role as transcriptional beta-catenin cofactors, fundamental for the transcription of Wnt target genes
Gene expression of muscle tissue in healthy controls
Bulk RNAseq was performed on muscle biopsies of 11 healthy matched controls, before and after, glucose and insulin infusion via hyperinsulinemic euglycemic clamps. In total we detected the expression of 20,849 in skeletal muscle. Principally, we see a clear distinction in gene expression patterns between basal and insulin states across all individuals (Figure 1). We observe 1,648 (837 upregulated and 811 down regulated) differentially expressed genes (p<0.01) between basal and insulin stimulated tissue biopsies in healthy controls (Figure 1b).
Notably highly significant up-regulated genes including CXCL1, IL6, IL1B, SLAMF1, and SOCS3 are cytokines with established functions in negatively regulating insulin signaling [25-28]. Highly downregulated include PDK4, KLF15, C1orf51, PIK3IP1, and TXNIP, the only gene shown to be differentially expressed in a similar study [29]. Insulin receptor substrate 1 (IRS1) and 2 were also found to be significantly downregulated.
Gene expression of muscle tissue in LRP5 mutants
Principally, we see a marred level of transcription in the individuals with LRP5 mutations. Individuals with LRP5 mutations demonstrated 916 differentially expressed genes (468 upregulated and 448 down regulated) (p<0.01) between basal and insulin stimulated tissue biopsies (Figure 1c). Seventy one percent of differentially expressed genes overlap with the differentially expressed genes in controls. Top upregulated genes in general had lower differential expression in comparison to controls, contrastingly, top downregulated genes were more differentially expressed (supplementary figure).
## Loading required package: ggrepel
Top 20 up-regulated genes across both mutant and control
Top 20 down-regulated genes across both mutant and control
Unique differentially expressed genes in LRP5
## 0.7055616 significantly differential genes overlap between mutant and wild typeHowever, many of these genes have a low logFC. Only 23 have an absolute change greater than 1, additionally gene enrichment reutrn no significantly enriched GO or REACTOME pathways, one KEGG pathway was returned relating to Gap junction
Gene enrichment
Gene enrichment was performed using gprofile. Gene were organised absolute logFC and submitted for enrichment via an ordered query allowing for genes with higher differential expression to have higher weight of enrichment.
KEGG and reactome contain more broad pathways and GO more detailed pathway. i.e One KEGG/REACTOME pathways may be multiple GO pathways.
GO analysis of the differentially expressed genes (Supplementary Data ?) demonstrate similarly enrichment in pathways relating to intracellular signaling, cell proliferation, vasculature, cellular, and structural development, although to a lesser degree (figure 1). However, in contrast to controls no significant enrichment in pathways regulating to insulin signaling are observed (figure 1). We identified genes within the mitogen-activated protein kinases (MAPKs), originally called extracellular signal-regulated kinases (ERKs) signaling pathway to be significantly differentially expressed between cases and controls (figure 1).
Further, MAPK is known to regulate the activities of several transcription factors, namely FOS and MYC both of which are greatly differentially expressed in controls in comparison to wild type.
Top KEGG pathways
Top REACTOME pathways
Overall, we see concordance between the top enriched KEGG and REACTOME pathways, both relating to interleukine and cytokine signaling, this is likely a reflection of negative feedback mechanisms switching insulin off. One interesting pathway is the resistance to EGFR tyrosine kinase inhibitors pathway. This is one of the few pathways more enriched in the mutant and effects three downstream pathways:
- rat sarcoma (RAS)/rapidly accelerated fibrosarcoma (RAF)/mitogen-activated protein kinase (MAPK) pathway
- phosphatidylinositol-3-kinase (PI3K)/protein kinase B (AKT) pathway and
- janus kinase (JAK)/signal transducers and activators of transcription (STAT) pathway, which stimulates mitosis, leading to cell proliferation and inhibition of apoptosis
How about the two main insulin signaling pathways? PI3K/AKT & MAPK/ERK
GO biological pathways functional enrichment
Overlapping Gene ontology terms were visualised and grouped using cytoscape.
Gene ontology (GO) analysis of the differentially expressed genes (Supplementary Data?) indicate highest enrichment in pathways relating to intracellular signaling, cell proliferation, vasculature, cellular, and structural development (figure 1). These pathways appear to be down stream of glucose uptake and insulin signaling within the muscle. Terms relating to insulin signaling were also significantly enriched, namely, response to peptide hormone (including insulin), regulation of protein phosphorylation and MAPK & ERK signaling pathways. The lower enrichment of insulin signaling related genes and pathways and the higher enrichment of other genes and pathways, most notably those involved in cytokine signaling which are an important negative feedback mechanism, suggests that the tissue is in the process of switching off signaling insulin signaling. This could be due to the timing of the muscle biopsies in relation to euglycaemic clamp stimulation suggesting that the biopsies were taken at the tail end of gene expression relating to insulin signaling and instead the tissue is shifting towards inhibition of insulin signaling.
Comparsions of pathways relating to insulin signaling
#### Genes enriched for peptide hormone response
#### Genes enriched for protein phosphorylation pathways
#### Genes enriched for JAK-STAT signalling
Muscle miRNA results
Bulk miRNA sequencing was performed on muscle biopsies of 12 controls and 10 individuals with LRP5 mutations. In total the expression of 358 miRNAs was measured. In contrast to gene expression we observe eight miRNAs were differentially expressed (P < 0.01) between cases and controls at basal state. Six of these same miRNA are significantly inversely differentially expressed in post euglycaemic clamp muscle tissue, indicating that these same miRNAs were degraded after glucose and insulin stimulation.
Contrastingly we see a large difference in miRNA expression is between wt and mut in the basal state, with 13 miRNAs (adjusted P < 0.05). These same miRNA show converse logFC in the mut basal vs insulin, indicating that these same mircoRNAs were used during insulin stimulation.
Degraded miRNAs
A caption
Association between M-value and differentially expressed miRNAs
Of these 8 differentially expressed miRNA how many are associated to M-values in the patients?
A linear model was used to test the association of normalised miRNA expression and M-values.
##
## Call:
## glm(formula = mir_counts[, "M_value"] ~ mir_counts[, "hsa-miR-224-5p"] +
## Age + Sex, family = "gaussian", data = mir_counts)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -112.84 -42.16 -17.87 53.92 136.67
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1130.580 197.127 5.735 3.94e-05 ***
## mir_counts[, "hsa-miR-224-5p"] -80.109 23.811 -3.364 0.004258 **
## Age -6.450 1.257 -5.133 0.000123 ***
## Sexmale 32.522 38.827 0.838 0.415393
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 5146.787)
##
## Null deviance: 231670 on 18 degrees of freedom
## Residual deviance: 77202 on 15 degrees of freedom
## AIC: 221.8
##
## Number of Fisher Scoring iterations: 2## [,1] [,2] [,3] [,4] [,5]
## (Intercept) 8.3618589 9.11242960 12.11958408 8.06648974 10.47398664
## mir_counts[, x] -0.3588533 -0.35282895 -0.56304684 -0.47288746 -0.40375099
## Age -0.0197690 -0.02300602 -0.02240703 -0.01326774 -0.01820745
## Sexmale 0.1350483 -0.02726950 0.14963232 0.10735380 0.14178982
## [,6] [,7] [,8]
## (Intercept) 8.88015254 4.37840364 7.8503496
## mir_counts[, x] -0.50740244 0.40885022 -0.2863919
## Age -0.02032587 -0.01430091 -0.0189073
## Sexmale 0.08595412 -0.01202168 0.1400002## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
Enrichment of miRNA target genes within differentially expressed genes
We then mapped these miRNA to gene targets and tested for enrichment within our differentially expressed genes, we found significant evidence of enrichment.
Gene targets of these miRNA are significantly enriched for differentially expressed genes identified.
## 3552 miRNA targets identified via seed pairing of 7 upregulated miRNAs## [1] "Results of enrichment test of miRNA targets between differentially expressed genes and randomly selected genes"##
## 2-sample test for equality of proportions with continuity correction
##
## data: c(sum(all_tar2$target %in% all_sig$Symbol), mean(rand_avg)) out of c(1856, 1856)
## X-squared = 62.709, df = 1, p-value = 2.396e-15
## alternative hypothesis: two.sided
## 95 percent confidence interval:
## 0.07666763 0.12760608
## sample estimates:
## prop 1 prop 2
## 0.2392241 0.1370873## [1] "Results of enrichment test of miRNA targets between control and LRP5"##
## 2-sample test for equality of proportions with continuity correction
##
## data: c(sum(all_tar2$target %in% sigWT$Symbol), sum(all_tar2$target %in% sigMUT$Symbol)) out of c(1651, 917)
## X-squared = 1.7016, df = 1, p-value = 0.09604
## alternative hypothesis: less
## 95 percent confidence interval:
## -1.000000000 0.006283277
## sample estimates:
## prop 1 prop 2
## 0.2356148 0.2595420We observed enrichment of targets to genes but are these genes being silenced as we would expected?
Lets compare the logFC of the targets versus non-targets for both wild and mutant
boxplot(sigWT$logFC_wt_InsVsBasal[sigWT$Symbol %in% all_tar2$target],
sigWT$logFC_wt_InsVsBasal[!sigWT$Symbol %in% all_tar2$target],
sigMUT$logFC_mut_InsVsBasal[sigMUT$Symbol %in% all_tar2$target],
sigMUT$logFC_mut_InsVsBasal[!sigMUT$Symbol %in% all_tar2$target],outline=F,
names = c("miRNA targets \n WT", "non-targets \n WT", "miRNA targets \n MUT", "non-targets \n MUT"),
title = "Comparsion of miRNA targets/non-targets LogFC")
text(x= 1.5, y= 1.5, labels= "0.02")
text(x= 3.5, y= 1.5, labels= "0.02")
miRNA conclusions
- We observe 8 differentially expressed miRNAs at basal level in LRP5 mutants
- These miRNAs are then degraded after euglycaemic clamp indicating a role in silecning genes
- We tested the association of these miRNAs to M-values and observe the miRNA (hsa-miR-224-5p) with the largest logFC also has the strongest inverse association to M-value. All other miRNAs although non-significant follow the expected direction based on expression
- Lastly, we tested if targets of these miRNA are enriched in our differentially expressed genes. We observe strong enrichment overall and suggestive enrichment for LRP5 differentially genes over wild type (p=0.09)